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ABSTRACT 


Magnetotellur ic  measurements  in  the  frequency  range 
0.001  to  1  Hz  were  made  at  the  University  of  Alberta 
Geophysical  Observatory  near  Leduc,  Alberta,  as  part  of 
the  crustal  conductivity  measurements  in  central  and 
southern  Alberta.  The  analog  data  obtained  were  converted 
into  digital  form  and  power  spectral  analyses  using  the  Fast 
Fourier  transform  techniques  were  made  on  the  digital  data. 
The  apparent  resistivity  curves  obtained  from  the  power 
spectral  density  estimates,  and  the  polarization  studies  of 
horizontal  electric  and  magnetic  fields  indicated  the 
presence  of  strong  resistivity  anisotropy  at  the  Observatory. 

Apparent  resistivity  curves  for  several  theoretical, 
anisotropic  layered  earth  models  were  computed  in  order  to 
match  with  the  experimental  curves.  The  results  thus 
obtained  indicated  that  the  resistivity  structure  beneath 
the  Observatory  corresponds  to  a  four  layered  anisotropic 
earth  model,  with  the  upper  isotropic  sedimentary  section 
underlain  by  two  anisotropic  layers  of  the  Precambrian  base¬ 
ment.  The  total  thickness  of  the  anisotropic  layers  and  the 
direction  of  the  major  axis  of  anisotropy  were  found  to  be 
10  km  and  N  50°  E  respectively.  A  resistivity  of  1000  ohm 
meters  was  obtained  for  the  bottom  layer. 
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CHAPTER  1 


INTRODUCTION 

During  the  summer  months  of  1966  magnetotellur. ic  measure¬ 
ments  were  carried  out  at  several  central  and  southern 
Alberta  stations,  in  order  to  delineate  the  conductivity 
structure  of  the  upper  crust  in  these  regions.  During  the 
same  period  continuous  registration  of  magnetotelluric  field 
variations  had  also  been  made  at  the  University  of  Alberta 
Geophysical  Observatory  (53°  13'  N;  113°  21'  W) ,  so  as  to 
facilitate  the  studies  on  magnetotelluric  source  fields.  The 
latter  measurements  were  also  intended  to  determine  the  con¬ 
ductivity  structure  at  the  Observatory,  which  is  to  be  the 
reference  station  for  future  magnetotelluric  studies  at  the 
University  of  Alberta.  The  subject  matter  of  this  thesis  is 
mainly  the  determination  and  interpretation  of  the  conductivity 
structure  at  the  Observatory,  from  the  magnetotelluric  data. 

1.1.  Geology  at  the  Measuring  Site 

The  Geophysical  Observatory  of  the  University  of  Alberta, 
where  the  present  magnetotelluric  measurements  we re  made,  is 
situated  in  the  central  plains  of  Alberta.  The  Precambrian 
shield  in  the  central  plains  is  overlain  by  about  2.5  km  of 

i 

flat-lying  sediments.  The  sedimentary  section  in  this  region 
has  been  fairly  well  studied  by  geophysical  and  geological 
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methods  in  the  course  of  exploration  for  oil.  A  typical 
cross-section  of  the  sediments,  with  well  log  resistivities 
is  illustrated  in  Figure  1. 

Garland  and  Burwash  (1959)  have  studied  the  petrology 
of  the  Precambrian  shield  using  well  log  samples  and  gravity 
anomalies.  Their  lithological  map  of  the  Precambrian  base¬ 
ment  is  depicted  in  Figure  2  on  which  is  superimposed  the 
location  of  the  Observatory.  Vozoff  et  al  (1963)  have 
carried  out  magnetotelluric  measurements  in  this  area 
and  found  that  no  significant  resistivity  anisotropy  exists 
in  the  sedimentary  section.  No  crustal  seismic  data  are 
available  for  this  location  but  the  seismic  crustal  studies 
by  Richards  and  Walker  (1959),  about  200  miles  farther  south, 
gave  depths  of  29  km  and  43  km  to  the  Conrad  and  Mohorovicic 
discontinuities.  They  also  predicted  heterogeneity  in  the 
upper  or  "granitic"  layer. 

1.2.  Historical  Review  of  the  Magnetotelluric  Method 

Almost  simultaneously  Tikhonov  (1950),  Kato  and  Kikuchi 
(1950)  and  Rikitake  (1950,  1951)  observed  the  existence  of  a 
correlation  between  the  geomagnetic  and  telluric  field 
variations  at  the  surface  of  the  earth,  and  investigated 
methods  for  exploring  to  great  depths.  However,  it  was  not 
until  Cagniard's  (1953)  classic  paper  that  the  method 
attracted  serious  attention.  Cagniard  presented  a  method  of 
interpreting  the  results  of  magnetotelluric  measurements  as 
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Elevation  A  B 


Figure  1:  Geological  cross-section  AB  (Figure  2) 
of  Edmonton  area,  with  typical  well-log 
resistivities  (after  Ellis,  1964) 
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Figure  2 : 


Lithology  of  the  Precainbrian  basement  as 
inferred  from  well  samples  and  gravity 
anomalies  (after  Garland  and  Burwash,  1959) 
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applied  to  geophysical  exploration.  Cagniard's  method 

4 

assumes  that  the  earth  is  horizontally  stratified  with  each 
layer  being  homogeneous  and  isotropic,  and  that  the  source 
fields  are  plane  electromagnetic  waves  impinging  on  the 
surface  of  the  earth.  Wait  (1954),  questioning  Cagniard's 
assumptions,  has  shown  that  Cagniard's  results  for  a  layered 
earth  are  valid  only  if  the  fields  themselves  do  not  vary 
appreciably  in  a  horizontal  distance  of  the  order  of  a  skin 
depth  in  the  earth.  Consequently,  the  field  should  be 
uniform  over  a  wide  area  to  permit  Cagniard's  interpretative 
procedure  to  be  applied.  Price  (1962)  has  indicated  that  the 
limitation  mentioned,  becomes  much  more  stringent  when  the 
magnetotelluric  method  is  applied  to  a  stratified  earth  in 
certain  important  instances.  However,  Madden  and  Nelson  (1964), 
considering  some  realistic  earth  models,  have  concluded  that 
the  plane-wave  assumption  is  valid  in  most  cases.  Srivastava 
(1965)  has  shown  that  the  effect  of  wavelength  of  the  source 
can  be  assumed  to  be  negligible  for  periods  less  than  1000  sec, 
when  determining  moderate  resistivity  distribution  in  shallow 
regions  (10  to  20  km)  within  the  earth. 

Following  Cagniard's  methods,  magnetotelluric  measure¬ 
ments  have  been  carried  out  successfully  by  several  investiga¬ 
tors  (Cantwell,  1960;  Wiese,  1962;  Srivastava  et  al .  ,  1963; 

Vozoff  et  al.,  1963;  Tikhonov  and  Berdichevskii,  1966; 

Hopkins  and  Smith,  1966;  and  others).  But  problems  have 
arisen  in  interpreting  the  magnetotelluric  data  in  areas  of 
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lateral  conductivity  variations.  Chetaev  (1960),  Kovtun 
(1961),  Rokityanski  (1961),  Cantwell  (1960),  Bostick  and 
Smith  (1962) ,  and  others  have  suggested  methods  to  interpret 
magnetotelluric  data  over  anisotropic  or  inhomogeneous 
earth.  Theoretical  studies  for  two  dimensional  inhomogeneities 
have  been  made  by  Neves  (1957),  d'Erceville  and  Kunetz  (1962), 
Rankin  (1962),  Weaver  (1963)  and  Mann  (1964).  Analog  model 
studies  over  two  and  three  dimensional  bodies  have  been  made 
by  Rankin  et  ai.  (1965),  Dosso  (1966)  and  others. 

1.3.  Outline  of  Thesis 

In  Chapter  2  electromagnetic  theory  applied  to  a 
layered  half-space  and  Cagniard's  magnetotelluric  effect  are 
studied.  The  magnetotelluric  theory  over  n-layer  anisotropic 
earth  models,  and  a  method  of  obtaining  apparent  resistivity 
curves  for  such  models  are  presented. 

A  brief  description  of  the  magnetotelluric  measuring 
equipment  and  the  techniques  of  handling  the  magnetotelluric 
data  are  given  in  Chapter  3. 

In  Chapter  4  polarization  diagrams  and  apparent  resis¬ 
tivity  curves  are  presented.  These  apparent  resistivity 
curves  are  interpreted  in  terms  of  a  four-layered  anisotropic 
earth.  An  interpretation  of  the  resistivity  structure  in 
terms  of  local  geologic  structures  is  attempted.  Chapter  4 
is  followed  by  conclusions  and  suggestions  for  further  work, 
references  and  three  Appendices. 
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CHAPTER  2 


MAGNETOTELLURIC  THEORY 

2.1.  Intrinsic  Impedance  of  a  Semi-infinite  Isotropic 
Half-Space 

For  the  case  of  plane  electromagnetic  waves  impinging 
on  a  homogeneous  isotropic  earth,  the  geometry  of  the  ionos¬ 
pheric  layers  and  the  frequencies  of  interest  (viz.  0.001  Hz 
to  1  Hz)  allow  the  effect  of  the  earth's  curvature  to  be 
neglected.  Because  of  the  enormous  contrast  between  the 
velocity  of  propagation  of  electromagnetic  waves  in  the  air  and 
in  the  earth,  the  refracted  wave  will  propagate  vertically 
downward  regardless  of  the  angle  of  incidence.  Furthermore 
the  component  of  current  across  the  surface  can  be  neglected. 

In  general  the  displacement  current  is  negligible  compared  to 
the  conduction  current  in  earth  materials,  in  the  frequency 
range  of  interest.  The  magnetic  permeability  in  most  cases 
of  interest  can  be  taken  to  be  that  of  free  space  (y  =  1) . 

If  in  addition,  one  considers  a  single  Fourier  component  of 
general  electromagnetic  disturbance,  the  time  dependence  can 
be  written  as  elwt  and  Maxwell's  equations  in  em  system  of 
units  can  be  written  as 

curl  H  =  4ttJ  =  4  ip  E 


(2.1) 


. 
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curl 

e 

=  -  i  a) 

“V 

H 

(2.2) 

div 

H 

=  0 

(2.3) 

and  the  Helmholtz 

(or 

diffusion) 

equation  for  the  electric 

•  "V  % 

field  E  m  the  medium  becomes 

2  -> 

,  2  -> 

(2.4) 

V  E 

k  E  =  0 

where  k^  = 

4 

7T  i  oj  a 

• 

(2.5) 

Considering  a 

cartesian  coordinate  system  with 

the 

x-y 

plane  on  the  earth 

1  s 

surface  and 

z  vertically  downwards, 

and 

assuming  that  the 

electric  vector 

E  in  the  plane  wave  is 

polarized  in  the 

X 

direction  and  the  magnetic  vector 

H  in 

the  y  direction, 

the  diffusion 

equation  (2.4)  can 

be 

written  as 

.2 

d  E 

X 

k  E  . 

(2.6) 

n  2 

dz 

X 

The  solution  of  Equation  (2.6)  in  the  case  of  a  semi-infinite 
medium  (z  -»•«>),  has  the  form 


(2.7) 


From  Equation  (2.2) 
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ltd 


dz 


(2.8) 
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Substituting  E  value  from  Equation  (2.7)  in  Equation  (2.8), 

A 

we  have 

Hy  =  TS  A  e'kZ  '  <2-9> 

Then  the  intrinsic  impedance  of  the  medium,  for  the  plane  wave, 
as  defined  by  Schelkunoff  (1943) ,  is 


Z 


E 

x 


(2  aT)  2 


-i  tt/4 
e 


(2.10) 


From  the  above  relation  one  can  see  that  the  intrinsic 
impedance  is  a  function  of  the  conductivity  of  the  medium  a, 
and  the  period  of  the  wave  T,  and  for  a  plane  wave  impinging  on 
a  homogeneous,  isotropic  semi-infinite  medium,  the  phase  of 
is  retarded  by  tt/4  radians  with  respect  to  that  of  E^. 

The  resistivity  of  the  medium  p  can  be  expressed  in  terms 
of  the  intrinsic  impedance  as 


(2.11) 


which  is  the  basic  equation  developed  by  Cagniard  (1953)  for 
the  magnetotelluric  method  of  exploration.  In  practical  em 
system  of  units  Equation  (2.11)  is  written  as 


P 


0.2  T 


(2.12) 


where 
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p 


resistivity  in  ohm-rue  ter  s 


E 


H 


T  =  period  of  magnetotelluric  disturbance  in  seconds 
x  -  disturbance  in  the  electric  field  in  iriv/km 
=  disturbance  in  the  magnetic  field  in  gammas 


The  depth  of  penetration  (or  skin  depth)  is  defined  as 
the  depth  p  in  kilometers,  where  the  amplitude  is  reduced  to 
fraction  1/e  of  what  it  is  on  the  surface  and  is  given  by 


(2.13) 


2.2.  Intrinsic  Impedance  of  an  n-layer  Homogeneous  Isotropic 


Half-space 


Cagniard  (1953)  obtained  formulae  for  the  intrinsic 
impedance  of  two  and  three  layered  models  and  provided  methods 
to  interpret  experimentally  observed  curves  in  terms  of  a 
two  layer  earth  model,  by"  curve  matching  techniques.  Yungul 
(1961)  presented  magnetotelluric  sounding  curves  over  three 
layer  earth  models  and  gave  interpretation  procedures. 
Srivastava  (1967)  published  a  catalog  of  apparent  resistivity 
and  phase  curves  over  two  and  three  layered  earth  models. 

Berdichevskii  (1960)  extended  Cagniard' s  theory,  for 


two  and  three  layered  earth  models,  to  an  n-layered  isotropic 
half-space  and  obtained  the  following  relation  for  the 
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intrinsic  impedance 


E 


Z  = 


x 


H 


z  =  0 


103 

k. 


coth  £  k^  +  arcth  [" 


where 


coth  (k2  + 


. .  +  arcth 


‘  kn-2 
L  kn-l 


coth  (k  n  +  arcth  — -r- —  ) 
n-1  n-1  k  . 

n 


..  ) 


(2.14) 


T_ 

k  .  =  wave  number  in  the  in  layer, 

m  1  ' 

P 

given  by  (-  4  it  i  u  /  p  )  2 

m 

h  =  thickness  of  the  m*'*1  layer 
m  J 

p  =  resistivity  of  the  mt^1  layer. 


It  becomes  somewhat  difficult  to  compute  intrinsic  impedance 
using  Equation  (2.14),  particularly  when  the  source 
dimensions  are  introduced  and  the  number  of  layers  in  the 
model  is  greater  than  three. 

In  order  to  compute  the  intrinsic  impedance  and  corres¬ 
ponding  apparent  resistivity  for  an  n-layered  isotropic  half¬ 
space,  avoiding  the  rather  laborious  methods  of  the  previously 
mentioned  authors,  a  straightforward  iterative  method  which 
produces  the  apparent  resistivity  curves  for  complex  models 
with  relative  ease  is  presented  here. 

Again,  considering  the  rectangular  coordinate  system 
with  the'  surface  of  the  earth  at  z  =  0  in  the  x-y  plane,  and 
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assuming  that  the  model  consists  of  n  homogeneous  isotropic 
plane  parallel  layers  underlain  by  a  homogeneous  isotropic 
half-space  (Fig.  3),  the  diffusion  Equation  (2.6).  in  each 
layer  is  written  as 


where 


,2 
a  E 


x,m 


? 

k~  E 

m  x  ,m 


(2.15) 


m 


4  TT  i  0) 


m 


(2.16) 


is  the  wave  number  for  the  m^  laver,  and  p  is  the  resis- 
tivity  of  the  m^  layer.  The  electric  and  magnetic  field  com¬ 
ponents  in  the  mt^1  layer  can  be  written  as  (Jones,  196  4) 

E  =  A  ekm(z-zm~l)  +  B  e“km  (z“zm--l^  (2.17) 

x,m  m  m 

where  A  and  B  are  constants.  By  using  Maxwell's  equations, 
mm  j  j  ir 

th 

the  horizontal  component  of  the  magnetic  field  in  the  m 

laver  H  can  be  written  as 
y,m 


-k 


H 


m 


y  ,m  ioj 


A  e 
m 


km(z  zm-l)  _  B  e“km(z"zm-l) 

m 


(2.18) 


Applying  the  boundary  condition,  viz.,  the  tangential  compo¬ 
nents  of  electric  and  magnetic  field  vectors  are  continuous 
at  the  boundary  z  =  zm,  one  gets  from  Equations  (2.17)  and 
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(2.18) 


k  h 

A  e  m  rn  +  B  e 
m  m 


k-mNn 


a  ^+1  hm 

Am+1  e  ^  Bm+1  e 


-k 


m+ 1  nm 


( ? 


k  h. 


•kmh, 


A  e“m“m  b  -  ‘  m 


m 


m 


“m+1 


m 


,  Km+1  hm  __  ^m+1  hm 

A  .  ,  e  B  ,  ,  e 

m+1  m+1 


.  (2 


Adding  Equation  (2.19)  to  Equation  (2.20),  one  gets 


*  +  ^m+1  „  ^m+1  ^m^Bm  ,  ^m  ^m+l 

Am  =  2  k  -  Am+1  e  +  - FT" 


m 


m 


D  -  (k  ,  +  k  )  h 

B  , ,  m+1  m  m 

m+1  e 


(2 


and  subtracting  Equation  (20)  from  Equation  (19) 
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^m  k"m+ 1  ( ^m+  ]_  + 

-Tk -  Am+1  e 

m 


^m  +  ^m+ 1  D  ^  ^m+ 1  Km  ^  ^m 

+ - 2 -  Bm+1  e 


(2 


If  the  lower  boundary  of  the  (n+1)  layer  extends  to 


(n+lf 


infinity  (hn+-^  ->  °°  )  ,  one  can  start  in  the  nether  world 


An+1  ”  °'  Bn+1 


1, 


(2 
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21) 
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and  we  have  very  simple  expressions  for  An  and  B  ,  viz. 


A 


n 


k  -  k  , 
n  n+1 


2  k 


n 


-  (k  ,  +  k  )  h 
n+1  n  n 


(2.24) 


B 


n 


k  +  k  , 
n  n+1 


2  k 


n 


-  (k  ,  -  k  )  h 
n+1  n  n 


(2.25) 


In  the  first  layer 


A. 


k2  +  kl 

2  k  ' 


*  ^(k2  kl)hl  ,  kl  k2  „  (k2+k1)h1 

A2  e  +  2  k,  B2  e 


(2.26) 


B. 


kl  k2 


2  k. 


(k2+ki )  hi  k2  +  kl  „  (k2  kl)hl 

A„  e  +  - y-Tj -  B0  e 

2  2  k-,  2 


(2.27) 


At  the  surface  (z  =  0) 


H 

y 


+  B- 


ioj 


(2.28) 

(2.29) 


The  intrinsic  impedance  for  an  n-layered  homogeneous  isotropic 
half-space  is  then  given  by 


Z  = 


iu) 


+  B 


i 


-  B 


1 


(2.30) 
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and  the  apparent  resistivity  in  practical  electromagnetic 
units  is 

pa  =  0.2  T 

Thus . starting  from  Equation  (2.23)  one  obtains  Equations 
(2.26)  and  (2.27)  by  utilizing  Equations  (2.21)  and  (2.22). 

Wait  (1954)  investigating  the  effect  of  non-plane  wave 
sources  on  magnetotelluric  theory,  has  shown  that  for  magneto- 

telluric  oscillations  of  frequency  0.1  Hz  on  a  homogeneous 

-  3 

earth  of  conductivity  a  =  10  mho/m,  the  spatial  variations 

of  E^  and  within  a  region  of  radius  approximately  35  km 

must  be  small,  otherwise  corrections  for  second  and  higher 

order  spatial  derivatives  of  the  source  fields  must  be  applied 

to  Cagniard's  plane  wave  theory  for  magnetotelluric  fields. 

He  has  also  pointed  out  that  if  the  micro-pulsation  sources 

were  in  a  lower  ionosphere  at  an  altitude  of  100  km,  corrections 

would  be  important  for  frequencies  less  than  0.1  Hz. 

Price  (1962)  has  investigated  the  effect  of  the  source 

dimensions  on  the  intrinsic  impedance  E  /H  .  For  a  homogeneous 

x  y 

earth  of  conductivity  a,  he  found  that 


E 


x 


H 


y 


(2.31) 


E 

X  ICO 


(2.32) 


where 


v  +  4  7i  i  a  a) 


(2.33) 
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and  v  is  a  separation  constant  which  has  a  dimension  L  '. 
According  to  Price  2tt/v  represents  the  linear 

dimension  of  the  source  field.  He  has  also  estimated  the 

-7  -1 

values  of  v,  and  found  that  they  range  from  1.6  x  10  cm 

-5  -1 

for  global  fields  to  1.6  x  10  cm  for  local  fields. 

Taking  into  account  the  dimensions  of  the  source  field, 
the  intrinsic  impedance  of  n-fayer,  homogeneous,  isotropic 
half-space  can  be  written  from  equations  (2.30),  (2.32)  and 

(2.33),  as 


Z 

s 


where  0  = 


ico 

0 


+  B 

-  B 


1 

1 


+ 


+  i 


+  v 


/ 


and 


4  TT  a  03  . 


(2.34) 


(2.35) 


Subscript  s  in  Equation  (2.34)  indicates  that  the  effect  of 
source  field  is  included.  The  apparent  resistivity  pg,  from 
Equation  (2.31)  and  (2.34)  is  written  as 

2 

pa  =  0.2  T  |  Z  |  (2.36) 

K  s  1  s  1 


A  program  has  been  written  for  the  IBM  360/67  digital 
computer  to  compute  apparent  resistivity  curves  for  n-layered 
earth  models.  One  such  model  together  with  apparent 


- 


. 


curve 
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Figure  4:  Computed  apparent  resistivity  curves  for  a  five  layered  earth 
model  for  various  values  of  source  dimension  v 
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resistivity  curves  for  various  values  of  v  is  depicted  in 
Figure  4,  It  may  be  noted  that  v  =  0  (the  source  is  infinitely 
large)  corresponds  to  Cagniard's  plane  wave  solution.  In  the 
rest  of  this  thesis,  v  is  taken  as  zero,  since  the  frequencies 
of  interest  and  depth  of  investigation  allow  the  effect  of 
source  dimension  to  be  neglected  (Madden  and  Nelson,  1964; 
Srivastava,  1965) . 


2.3.  Intrinsic  Impedance  of  Anisotropic  Media 

The  scalar  impedance,  as  developed  in  the  previous 
section,  allow  interpretations  only  in  terms  of  homogeneous 
and  isotropic  earth  models.  Cantwell  (1960)  developed  a 
procedure  for  determining  the  effects  of  anisotropy  and/or 
two  dimensional  inhomogeneous  structures  in  the  earth  in  terms 
of  an  impedance  tensor  of  the  form 


E 

X 

E 

y 

H 

X 

H 

y 

(3.27) 


The  electric  field  in  one  direction  may  depend  on 
magnetic  field  variations  parallel  to,  as  well  as  perpen¬ 
dicular  to,  its  direction.  A  method  of  calculating  the 
tensor  elements  involves  computing  the  Fourier  components  of 
E  and  H  for  two  independent  observations,  then  solving 
Equation  (2.37)  for  both  the  observations  simultaneously  for 
the  four  elements.  Swift  (1967)  has  made  an  extensive  study 


Jo 
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on  the  characteristics  of  theoretical  and  measured  impedance 
tensors . 

Bostick  and  Smith  (1962)  have  used  the  admittance  tensor, 


and  found  that  a  rotation  of  the  measuring  axis  to  minimize 
the  main  diagonal  terms  in  the  impedance  tensor  yields  a 
reasonably  consistent  estimate  of  the  direction  and  magnitude 
of  a  major  inhomogeneity  in  the  apparent  resistivity  profile. 

Mann  (1965)  has  analyzed  the  problem  of  plane  wave 
incidence  on  an  anisotropic  earth,  and  pointed  out  there  exist 
two  skin  depths  (Appendix  1)  corresponding  to  two  principal 
directions  of  anisotropy,  and  the  impedance  functions 
obtained  at  the  surface  are  highly  dependent  on  polarization 
as  well  as  on  anisotropy. 

O'Brien  and  Morrison  (1967)  have  studied  the  behaviour 
of  electromagnetic  fields  in  an  n-layered  anisotropic  half¬ 
space,  and  have  obtained  the  following  equations  for  the 

t  h 

total  components  of  electric  and  magnetic  fields  in  the  m 
layer,  where  the  axes  of  anisotropy  lie  in  the  x-y  plane. 

The  geometry  is  shown  in  Figure  5. 
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Figure  5 


Relationship  of  conductivity  anisotropy 
to  coordinate  axes 
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(2.39) 


and  the  plus  or  minus  superscripts  indicate  propagation  in 
the  plus  or  minus  z  direction  respectively,  where  subscripts 
1  and  2  refer  to  the  modes  of  propagation  corresponding  to  the 
two  principal  directions  of  anisotropy. 

From  Equation  (2.38)  the  boundary  conditions  may  be 
expressed  in  matrix  form,  as 
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and  T  is  identical  except  that  all  superscripts  are 

changed  to  m-1  and  the  subscripts  on  6  only  are  changed  to 

m-1 . 

The  electric  field  components  in  the  m-lth  layer  may 
be  determined  from  Equation  (2.41) ,  and  they  may  be  written 
in  matrix  form,  as 


where 
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Am  , 

3  1 1 
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(2.42) 


(2.43) 


For  the  (n+1)  layer  model  of  Figure  3,  the 
electric  fields  are  given  in  terms  of  the 
into  the  half-space,  by 
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and  A  . 

1  3  f  1 


is  defined  in  Equation  (2.43). 


The  surface  magnetic  fields  can  be  computed  from  the  surface 
electric  fields  by  Equations  (2.38).  Then  the  impedances 
and  hence  the  apparent  resistivities  in  the  x  and  y  directions 
can  be  calculated  using  Equations  (2.10)  and  (2.12)  respectively. 

Due  to  the  computational  difficulties  involved  in  a 
general  solution,  O'Brien  and  Morrison  (1967)  considered  only 
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the  case  where  0  is  a  constant,  independent  of  depth. 

Since  the  two  modes  of  propagation  for  a  single  layer 
anisotropy  are  independent  in  the  directions  of  principal 
axes,  it  is  valid  to  treat  a  single  layer  anisotropy  or  a 
multilayered  anisotropy  for  which  the  principal  directions 
are  the  same,  in  a  straight  forward  manner  with  computational 
ease,  to  obtain  the  impedances,  and  hence  the  apparent 
resistivities  at  the  surface. 

Considering  a  model  consisting  of  an  isotropic  upper 
layer,  two  intermediate  anisotropic  layers  and  isotropic 
bottom  layer  extending  to  infinity  (Figure  6)  and  representing 
the  resistivities  along  the  principal  axes  1  and  2  in  each 
layer  by  p ^  and  p2,  the  electric  and  magnetic  field  components 
and  E^,  and  H 2  can  be  obtained  by  using  the  formulae 

3.  3 

(2.28)  and  (2.29).  The  apparent  resistivities  p^2  and  p21  in 
the  two  principal  directions,  are  then  given  by  Equation 

3  3 

(2.31).  The  apparent  resistivities  p  and  p  along  any 
measuring  axes  x  and  y  may  then  be  obtained  by  the 
relation  (2.31) 
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Figure  7:  Relationship  of  electric  and  magnetic 
field  components  at  the  surface  of 
anisotropic  earth  (1  and  2  are  the 
principal  directions  of  anisotropy) 


(2.48) 
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and  0  is  the  angle  to  the  axis  of  measurement  x  in 
clockwise  direction  from  the  major  axis  of  anisotropy 
The  geometry  is  shown  in  Figure  7  and  a  set  of  master 


the 

1. 

curves 


obtained  for  the  model  considered  here  is  given  in  Figure  8. 
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CHAPTER  3 

DETECTION  AND  ANALYSIS  OF  MAGNETOTELLURIC  SIGNALS 
3.1.  Detection 

The  geographical  north-south  and  east-west  components 
of  variations  in  the  geomagnetic  and  telluric  fields  are 
recorded.  A  schematic  representation  of  the  magnetotelluric 
recording  system  is  shown  in  Figure  9. 

The  telluric  field  variations  are  derived  from  the 
potential  differences  between  copper  electrodes  spaced 
approximately  500  ft  (north-south)  and  340  ft  (east-west) 
apart.  A  100  y  f  capacitor  is  introduced  into  the  circuit  to 
block  DC  potentials  due  to  polarization  at  the  electrodes. 

This  capacitor  and  the  total  input  resistance  to  the  opera¬ 
tional  amplifier  summing  junction  produces  the  desired  lead 
corner  at  0.0073  Hx  (Figure  10). 

Induction  coils  consisting  of  32000  turns  of  wire,  wound 
on  a  one  inch  diameter,  5  feet  long,  high  permeability  core 
are  used  to  measure  variations  in  the  geomagnetic  field.  The 
calibrated  responses  of  electric  and  magnetic  systems  are 
shown  in  Figure  10. 

The  recordings  have  been  obtained  in  analog  form  on 
seven  channel  FM  magnetic  tapes  using  Precision  Instruments 
tape  recorder,  at  15/16  i.ps .  Simultaneously  six  channel 
brush  paper  recording  is  also  carried  out  in  order  to  provide 
an  editing  facility  prior  to  digitizing.  Time  signals  from 
W.W.V.B.  are  recorded  on  the  seventh  channel. 
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Block  diagram  of  magnetotelluric  recording  system 
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3.2.  Digital  Processing 

Records  approximately  three  hours  in  length  and  free 
from  discontinuities  are  selected  for  digital  conversion. 

A  block  diagram  of  the  Analog  to  Digital  Conversion  system 
is  given  in  Figure  11.  The  signal  in  analog  form  is  alias 
filtered  and  sampled  by  the  Analog  to  Digital  Converter 
at  a  rate  of  15  samples  per  second  per  channel.  This  gives 
a  digitizing  interval  of  32/15  seconds,  with  Nyquist  frequency 
of  approximately  0.25  Hz. 

Digitized  data  in  11  bit  binary  form  are  recorded  in  blocks 
of  1200  two  character  words  on  seven  channel  asynchronous 
Kennedy  digital  tape  recorders.  Blocks  are  read  alternately 
in  each  of  the  two  output  tape  recorders  (Figure  11)  in  order 
to  provide  the  necessary  Inter  Record  Gap  required  for  input 
to  the  IBM  360/67  digital  computer. 

3.3.  Spectral  Analysis 

In  applying  the  magnetotelluric  method  to  compute  the 
apparent  resistivity  of  the  earth  from  Equation  (2.12),  the 
ratio  Ex/Hy  or  Ey/Hx  as  a  funct:i-on  of  period  (or  frequency) 
can  be  obtained  by  visual  inspection  of  a  number  of  clear 
sinusoids  (Berdchevskii ,  1960;  Srivastava  et  al.,  1963)  or 
from  power  spectral  density  estimates  (Cantwell,  1960; 

Ellis,  1964;  Hopkins  and  Smith,  1966;  and  others).  The 
latter  method  is  adopted  in  the  present  studies. 
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Figure  11:  Block  diagram  of  Analog  to  Digital  Conversion  system 
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Cantwell  (1960)  has  considered  the  magnetotelluric  field 
components  as  samples  of  a  stationary  stochastic  process  and 
applied  the  numerical  methods  of  Blackman  and  Tukey  (1958) 
to  the  data.  He  has  also  shown  that  Equation  (2.12)  of 

Cagniard  may  be  expressed  in  terms  of  power  density  spectra, 
as 

a 
P 

where 

a 
P 

T 

P 

E  (f) 
x 


A  brief  description  of  the  spectral  computational  technique 
employed  in  the  present  studies  is  given  below. 

The  spectral  analysis  is  based  on  the  assumption  that 
magnetotelluric  records  are  from  stationary  processes  with 
zero  mean  values.  If  the  mean  value  is  not  zero,  then  the 
power  spectral  density  function  will  exhibit  a  large  peak  at 
zero  frequency  which  will  distort  estimates  at  other  frequencies. 
For  this  reason  the  original  records  are  conditioned  by  sub¬ 
tracting  the  sample  mean  values.  A  second  correction  is 


0 . 2  T 


E  (f ) 
x 

H  (f) 

y 


(3.1) 


apparent  resistivity  in  ohm  meters 
period  in  seconds 

power  spectral  density  of  telluric  field 

.  2 

component  E  in  (mv/km)  Hz 

X 

power  spectral  density  of  the  geomagnetic 

2 

field  component  H  in  (gammas)  /IIz. 
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applied  to  correct  for  the  slowly  varying  linear  trend  (that 
is,  non  zero  slope  with  respect  to  time).  This  may  be  due 
to  instrumentation  drift  or  to  an  actual  underlying  linear 
trend  in  the  data,  such  as  diurnal  variations. 

If  x(t)  is  the  original  series  of  length  R^r  the 
corrected  series  x(t)  is  written  according  to  Bendat  and 
Piersol  (1966) ,  as 

R1 

x ( t )  =  x ( t )  -  x  -  a  (t  -  ~y)  (3.2) 

ILL  A  Z 


0  <  t  Rn 
^  1 


where 


x . 
i 


(3.3) 


denotes  the  sample  mean  value  of  x(t)  over  (0,  R-^ )  ,  the 

parameter  a  denotes  the  average  slope, 

x 


a 


x 


1 

v(N- v  ) 


N 

n  =  N 


(3.4) 


where  N  is  the  total  number  of  data  samples  and  v  is  the 
largest  integer  less  than  or  equal  to  N/3. 

The  Fast  Fourier  Transform  (FFT)  algorithm  of  Cooley 
and  Tukey  (1965)  (Appendix  2)  is  used  to  compute  the 
covariance  functions  and  spectral  density  estimates. 
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The  Discrete  Fourier  Transform  (DFT )  of  a  series  x(t) 

4 

is  defined  as 


X  (W) 


N  -  1  ~2  l  1  t  W 

SI  x ( t )  e  N 

t  =  0 


W=0,  1,  2...N-1  (3.5) 


where  W  is  the  digital  frequency  index  and  N  is  the  total 
.number  of  data  points. 

Covariance  functions  are  obtained  employing  the  convolu¬ 
tion  theorem  techniques  of  Sande  (1965)  rather  than  summing 
the  lagged  products  which  is  computationally  expensive.  For 
example,  the  autocovariance  function  of  the  series  x(t)  is 
given  by 


N '  - 1 

1  X*  (W)  X(W) 

*  “2  7T  i  W  T 

K  e  N' 

o 

II 

L  J 

where  t  is  the  lag  number,  N'  is  the  total  number  of  points 
including  zeros  added  to  the  series  x(t)  to  maintain  the 
cyclic  nature  of  the  data,  X (W)  is  the  DFT  of  x (t)  as 
defined  in  Equation  (3.5)  and  *  represents  the  complex 
conjugate . 

The  raw  covariance  functions  thus  obtained  are  smoothed 
using  the  lag  window  of  Parzen  (1961)  with  the  weighing 


coefficients 


■ 
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W(,)  =  1-6  (I)2  +  6  (I)3 


I 


<  M 


=  0 


t  >  M 


(3.7) 


where  M  is  the  maximum  number  of  lags.  The  Parzen  lag 
.  window  has  proved  to  be  more  efficient  for  magnetotelluric 
analysis.  It  involves  less  computational  effort  than,  for 
example,  the  Daniel  window,  at  the  same  time  giving  positive 
power  estimates  with  extremely  low  side  lobes,  and  further, 
it  keeps  the  sample  coherence  between  its  theoretical  limits 
±  1. 

The  power  density  spectra  are  obtained  by  taking  the 
Fourier  transforms  of  the  smoother  covariance  functions.  The 
autopower  density  spectrum  of  the  x(t)  series  is  thus  given  by 

M  =  1  -2  tt  i  W  t 

P  (W)  =  yZ.  a  (t)  e'  M  (3.8) 

x v  '  n  m 

T  =  U 


where  a  (t)  is  the  smoothed  (or  modified)  value  of  a(x)  in 
m 

Equation  (3.6) . 

The  crosscovariance  function  of  x(t)  and  y(t)  series  is 
defined  as  the  Fourier  transform  of  the  product  X* (W)  Y (W) , 


1— I 

1 

l _ 

r 

+  -2  TT  i  w  T 

n 

X*  (W)  Y(W) 

e  N ' 

W  =  0 
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- 

*  .  »  •  V  , 
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where  ±  (W)  is  the  DFT  of  the  y (t)  series  and  the  rest  of  the 
nouation  is  the  same  as  in  Equation  (3.6) .  The  crosspowsr 
density  spectrum  is  then  obtained  by  taking  the  Fourier 
transform  of  the  smoothed  crosscovariance  function. 

The  sample  coherence  between  orthogonal  electric  and 

v 

magnetic  field  components  is  also  computed  using  the  definition 


Y  (f ) 


where 


P  (f ) 
xy  


/VfTpy(f) 


(3.10) 


Y  (f ) 


sample  coherence 


P  (f ) 
xy 


P  (f),P  (f) 
x  '  y 


crosspower  density  spectrum 


autopower  density  spectra  of  x  and  y  series 


Examples  of  the  autopower  density  spectra  of  orthogonal 
electric  and  magnetic  field  components  and  the  sample  coherence 
between  them  are  given  in  Figures  12  and  13.  The  number  of 
points  and  the  maximum  number  lag  used  were  4096  and  512 
respectively . 

If  the  spectral  density  estimates  at  a  frequency  f  are 

to  be  useful,  only  frequencies  close  to  f  should  have 

substantial  contribution.  If  the  spectrum  rises  steeply  as 

f  leaves  f  ,  an  extremely  rapid  fall-off  in  the  spectral 
o 

window  is  required.  To  avoid  this  problem  Blackman  and  Tukev 
(1958) ,  Ellis  (1964)  and  others  have  recommended  the  pre¬ 
whitening  of  the  data  before  computing  the  spectrum. 
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Figure  12:  Autopower  density  spectra  of^_Ex  (north 

south  electric)  and  Hy  (east-west 
magnetic)  and  the  sample  conerence 
between  them 
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Figure  13:  Autopower  density  spectra  of  Ey 

(east-west  electric) and  Hx  (north- 
south  magnetic)  and  the  sample 
coherence  between  them 
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In  order  to  test  the  whiteness  of  the  present  data,  a 
prewhitening  filter  of  the  form 

zi  =  x..  +  K  x±  _  ±  (3.11) 

with  K  =  -0.9  is  applied  to  original  data  before  computing  the 
spectrum  and  the  calculated  spectrum  is  then  divided  by  the 
power  transfer  function, 

Y  (f )  =  (1  +  K2  +  2K  cos  2  tt  fAt)  ,  (3.12) 

of  the  prewhitening  process.  It  was  found  that  this  prewhiten¬ 
ing  operation  did  not  significantly  improve  the  spectral 
estimates  for  the  range  of  frequencies  considered  in  this 
thesis.  This  may  be  due  to  the  fact  that  a  considerable  amount 
of  prewhitening  is  achieved  in  the  design  of  the  recording  system 
and  further,  the  spectra  were  relatively  flat. 
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CHAPTER  4 

4 

RESULTS  AND  INTERPRETATION 

•  •  •  3.  3 

Apparent  resistivity  curves  p  and  p  in  the  north- 

2  xy  yx 

south  and  east-west  directions  respectively,  were  obtained 
following  the  power  spectrum  analysis  methods  outlined  in 
Chapter  3,  and  using  Equation  (3.1).  Spectral  density 
estimates,  for  which  the  sample  coherence  between  orthogonal 
electric  and  magnetic  field  components  was  greater  than  0.9, 

3  3 

were  used  to  compute  apparent  resistivities  p  and  p 

L  x  ^  xy  yx 

The  average  apparent  resistivity  values  thus  obtained  from 
20  "good"  data  sets,  recorded  on  different  days  and  times  at 
the  University  of  Alberta  Geophysical  Observatory,  are 
plotted  as  functions  of  period  in  Figure  14.  The  error  bars 
in  the  figure  correspond  to  maximum  deviation  from  the  mean. 
These  apparent  resistivity  curves  deviate  significantly  for 
periods  greater  than  20  seconds,  suggesting  a  strong  anisotropy 
at  the  measuring  site.  Similar  types  of  anisotropy  are 
reported  by  Fournier  (1962)  at  Garchy,  Nivre  in  France,  and 
by  Hopkins  and  Smith  (1966)  at  Pecos,  Texas. 

Polarization  studies  have  also  been  made  ^for  all  data 
sets  chosen  for  the  analysis.  Plots  of  the  electric  and 
magnetic  field  polarization  values  are  line  printed  directly 
from  the  computer  and  are  illustrated  in  Figures  15  and  16 . 
These  plots  which  are  typical  of  all  the  data  sets  are  each 
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calculated  on  the  basis  of  3600  digital  values.  These 
polarization  diagrams  are  fitted  to  ellipses  and  the  parameters 
of  the  ellipses  viz.,  the  azimuth  of  the  major  axis  with 
respect  to  the  geographical  north  (X-axis)  and  the  ellipticity 
are  obtained  by  using  the  following  equations 

n 


a 


1 

2 


tan 


xj.  yi 


i=l 


n 

2: 

i=l 


n  2  n 

1-  ( ?-  xi  /  i ) 

'i=l  i=l  1; 


2  .  2 
y_.  x  tan  a 


(4.1) 


(4.2) 


(  21  x .  /  y .  )  -  tan  a 

i=l  1  i=l  1 

where  a  is  the  azimuth  of  the  major  axis  of  the  ellipse 
measured  in  clockwise  direction  with  respect  to  North,  and  a 
and  b  represent  the  major  and  minor  axes  of  the  ellipse 
respectively.  The  parameters  of  the  ellipses  thus  obtained 
are  given  in  the  corresponding  figures.  The  non-orthogonality 
of  the  resultant  horizontal  electric  and  magnetic  fields  for 
the  Pc4  type  pulsations,  which  dominate  the  results,  is  further 
evidence  of  anisotropy  at  this  location. 

Various  anisotropic  earth  models  (Appendix  3)  were  considered 
according  to  the  theory  outlined  in  Chapter  2,  in  order  to  get 
a  resonablv  good  fit  with  the  observed  apparent  resistivity 


curves . 


A  model  which  yields  apparent  resistivity  curves 
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which  are  in  clcse  agreement  with  observations,  is  shown  in 

Figure  17.  A  set  of  resistivity  curves  obtained  for  the 

model  given  in  Figure  17,  is  shown  in  Figure  18  for  various 

values  of  0,  the  angle  from  the  major  axis  of  anisotropy  in 

clockwise  direction  to  the  north-south  axis  of  measurement. 

The  p  and  pa  curves  for  0  =  -50°  give  a  reasonable  fit  to 

yx 

the  observed  curves  as  illustrated  in  Figure  14. 

The  period  20  seconds,  at  which  the  observed  pa  and  pa 

xy  yx 

curves  start  to  deviate  significantly  from  one  another, 
corresponds  to  a  depth  of  approximately  6  km.  Thus  it  is 
reasonable  to  assume  that  the  sedimentary  section  which  is 
approximately  2.5  km  thick  is  an  isotropic  single  layer  for 
modelling  purposes.  An  average  resistivity  of  8  ohm-meters 
is  assigned  to  this  layer  on  the  basis  of  well  logs  (Figure  1) 
and  previous  magnetotelluric  measurements.  This  layer  corres¬ 
ponds  to  the  first  layer  of  the  model.  The  upper  part  of  the 
basement  is  mainly  gneissic  in  composition  at  the  measuring 
site,  however,  there  are  significant  variations  in  lithology 
over  short  distances  (Figure  2).  The  parameters  viz . , 
resistivities  and  thicknesses  for  the  second,  third  and 
fourth  layers  of  the  model  are  selected  in  order  to  fit  the 
experimental  curves.  The  low  resistivity  values  chosen  for 
the  third  layer  do  not  agree  with  the  results  of  the  previous 
workers,  except  for  the  result  by  Vozoff  et  al.  (1963)  at 
Kavanagh  for  which  they  did  not  attempt  an  interpretation. 

The  apparent  resistivity  curves  obtained  by  them  for  the 
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Four  layered  anisotropic  earth 
model 
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Kavanagh  station,  which  is  a  few  miles  south  of  Leduc,  are 
shown  in  Figure  19.  The  depth  to  which  the  lithological  units, 
shown  in  Figure  2,  extended  within  the  basement  was  given  by 
Garland  and  Burwash  (1959)  from  gravity  studies.  The  depth 
quoted  by  them  w as  8  km,  which  is  in  good  agreement  with  the 
total  thickness  of  the  anisotropic  layers  in  the  model  presented 
here.  The  1000  ohm  meters  resistivity  of  the  bottom  layer  can 
be  varied  considerably  with  little  change  in  the  shape  of  the 
apparent  resistivity  curves,  however,  only  an  extension  of 
measurements  at  the  long  period  end  of  the  spectrum  could 
provide  conclusive  results. 

The  major  direction  of  anisotropy  which  is  N  50°  E  cannot 
be  interpreted  unambiguously  with  the  present  state  of  know¬ 
ledge  about  the  basement  structure.  This  anisotropy  may  be 
related  to  major  geologic  structures  such  as  the  Rocky  Mountains 
or  to  a  more  local  conductivity  anomaly.  The  effect  of  distant 
vertical  discontinuities  on  measured  surface  apparent 
resistivity  has  been  shown  by  Rankin  et  al.  (1965)  from  their 
analog  model  studies. 

The  direction  of  the  electric  field  vector  which  is 
predominantly  polarized  N  58°  W  is  close  to  the  value  obtained 
by  Srivastava  et  al.  (1963)  at  Cooking  Lake  which  is  a  few 
miles  north-east  of  the  Observatory.  They  interpreted  this 
north-west  direction  of  the  electric  vector  as  due  to  the 
structural  trends  in  the  sedimentary  layers.  An  examination 
of  analog  records,  however  indicates  that  the  contribution  to 
the  electric  field  polarization  (Figure  15)  comes  mainly  from 
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the  periods  ranging  from  60  to  150  seconds,  which  have 
much  deeper  penetration. 
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CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  WORK 

From  the  apparent  resistivity  curves  presented  in  this 
thesis,  it  may  be  seen  that  employing  refined  instrumentation 
and  "good"  statistics  reduces  the  scatter  in  the  apparent 
resistivity  curves  to  a  level  such  that  reliable  comparison 
can  be  made  with  model  calculations. 

While  neither  crustal  seismic  data  nor  aeromagnetic 
data  are  available  for  this  particular  location,  a  pronounced 
directional  trend  is  noticeable  along  the  direction  of  the 
major  axis  of  anisotropy  in  the  aeromagnetic  map*  of  the 
Cooking  Lake  area,  which  is  a  few  miles  North  of  the 
Observatory.  In  order  to  make  a  satisfactory  interpretation 
of  the  apparent  anisotropy  observed  at  the  Observatory  a 
close  spaced  magnetotelluric  network  with  other  geophysical 
coverage  is  necessary.  Further,  an  analysis  of  the  vertical 
component  of  the  magnetic  field  and  the  determination  of 
phase  relationships  between  horizontal  electric  and  magnetic 
field  components  may  be  helpful  in  interpreting  this 
resistivity  anisotropy. 

Care  must  be  taken  to  avoid  the  problem  of  measuring  field 
components  in  the  direction  of  polarization.  This  can  be 
done  bv  reconnaissance  survey  in  order  to  determine  the  appro¬ 
priate  direction  for  the  sensors  or  by  recording  two 
orthogonal  components  at  various  orientations. 


*  Geophysics  paper  33,  G.S.C.,  OTTAWA,  1951. 
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APPENDIX  1 


Al.l.  Anisotropic  Half-Space 

For  a  half-space  whose  conductivity  is  different  in 
different  directions,  Ohm's  law  is  written  as 

|  3  |  =  |  <t  l  I  £  |  .  (Ai.i) 

1  1  1  m,  n  1  1  1 

For  a  plane  wave  incidence 

J  =  0  (A1.2) 

z 

and  following  the  cartesian  notation  used  in  Figure  5 
(Chapter  2) 

J  =  an  E  ,  cos  0  +  a0  E  .  sin  0 

x  lx'  2  y 

J  =  -a-,  E  .  sin  6  +  a0  E  .  cos  0.  (Al.3) 

y  lx  2  y 

The  E  ,  and  E  .  can  be  written  in  terms  of  E  and  Ey 
x '  y  J 

E  ,  =  E  cos  0  -  E  sin  0 

x '  x  y 

E  ,  =  E  sin  0  +  E  cos  0  (Al.4) 

y '  x  y 


From  Equations  (Al.3)  and  (Al.4) ,  we  write 

j  =  a,  E  cos2e  -  a,  E  cos0  sin0  -  E  sin20  +  a2  E  sin9  cos0 
x  i  x  1  y  *  7 

2  2 
j  =  -<j  E  cos 0  sin0  +  o,  E  sin  0  +  a  E  sin0  cos0  +  a  E  cos  ( 
u  v  i  x  1  y  ^  ^  j 


(Al. 5) 


Ill 


or,  in  matrix  form 


3 


(A1.6) 


where 
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2  .  2 

cos  0  +  o' 2  sin  0  ^a2~°]  ^  s:*-n0  cos0 


(°2-ai^  sin0  cos© 


.  2  2 
sm  0  +  cos  0 


(A1.7) 


Neglecting  the  displacement  currents  and  taking  the  magnetic 
permeability  of  the  medium  as  unity,  from  Maxwell's  equations 
(in  em  system  of  units) 


k 

E 

y 
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e 
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E 
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CO  H 

y 

-i  k 
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y 

4  71  (ai]_ 

E 

X 

+ 

an  0  E  ) 

12  y 

i  k 

H 

X 

4  71  (o12 

E 

X 

+ 

o  22  Ev)  (Al.8) 

where  i  : 

co  is  the 

radian 

frequency  and  k  is  the 

propagation  constant. 

Solving  for  k  from  the  set  of  equations  (Al.8),  we  have 


IV 


t,2  _  4  it  i  a) 

(on 

+  a22) 

(011  + 

a22 )  +  4  (a12a21 

K  2 

"T 

“  all 

a22 ' 

h 

(A1.9) 

Substituting 

al,  2  s 

from 

(Al. 

7 )  ,  we 

get 

ki  = 

±  (4 

x  h 

it  a  ^  a) ) 

0i  ( tt/4  ) 

k2 

±  (4 

x  % 

7T  0  2  oj  ; 

e1  ( m / 4 )  (A1.10) 

Thus  there  are  two  independent  modes  of  propagation 
corresponding  to  two  principal  directions  of  anisotropy.  Then 
the  two  skin  depths  are  given  by 

p^  =  (2  n  u)  2 

P2  =  (2  it  c>2  w)  2 


(Al.ll) 


»  T 


' 


* 


V 


APPENDIX  2 

A2.1.  Fast  Fourier  Transform  Algorithm 

The  Fast  Fourier  transform  (FFT)  algorithm  is  a  method 
of  computing  the  finite  Fourier  transform  of  a  series  of  N 
(complex)  data  points  in  approximately  N  log2  N  operations. 

The  history  (Cooley  et  al.  1967)  of  this  algorithm  can  be 
traced  back  to  the  paper  of  Danielson  and  Lanczos  in  1942. 

But  it  was  not  until  Cooley  and  Tukey's  (1965)  paper  it  has 
attracted  any  attention.  Following  Cooley  and  Tukey's  paper, 
Gentleman  and  Sande  (1966),  Bingham  et  al.  (1967)  and  Welch 
(1967)  have  pointed  out  the  usefulness  of  this  algorithm  for 
a  variety  of  problems,  such  as  the  calculation  of  auto  and 
crosscovariances,  auto  and  crosspcwer  spectra  etc.  The 
basic  algebra  of  the  Fast  Fourier  transform  is  described  here. 

The  Discrete  Fourier  transform  (DFT )  of  a  time  series 
x(t)  is  defined  as 

(-2  tt  i  t  W) 

X(W)  =  ~T0  x(t)  e  N  (A2.1) 

where  N  is  the  number  of  points. 

For  algebraic  convenience  we  write 

e2  ”  i  Y  =  e(Y)  .  (A2.2) 


e(Y)  =  1  if  Y  is  an  integer. 


VI 


Equation  (A2.1)  can  be  written  as 
N-l 

X(W)  =  21  x(t)  e  (~) 

t=0  N 


(A2.3) 


Supposing  N  has  factors  G  and  H  so  that  N  =  GH,  and 
writing  W  =  a’  +  b'G  and  t  =  b  +  aH,  Equation  (A2.3)  can  be 
transformed  to 


X (a ' +b ' G)  = 


H-  1  G- 1 

YZ  x  (b+aH)  e  (" - — !°-+-aK — ) 


b=0  a=0 


GH 


IT  IT  x  (b+aH)  e('^|  -  21®  -  L'b  -  ab,) 

b=0  a=0  Gh  G  H 


H-l  G-l 


=  E  r  *  (b+aH)  eC^)  ef^')  e  (-2§-l) 


.--aa 


-bb 


b=0  a=0 


GH 


H 


and  since  ab 1  is  an  integer  so  that  e(ab')  =  1 


H- 1  ,  ,  ,  , ,  G- 1  . 

YZ  e(  -77")  e  (  — — )  2l__  X  (b+aH)  e(  -^-) 

b=0  b  a=0  G 


(A2.4) 


To  obtain  the  Fourier  transform  of  the  complete  series,  the 
inner  sum  of  length  G  need  be  evaluated  for  only  H  different 
values  of  b,  and  G  different  values  of  a'  requiring  G.HG 
multiplications;  then  the  outer  sum,  which  is  of  length  H 


• 

VI 1 


need  be  computed  for  H  values  of  b',  and  G  values  of  a'  using 

H.GH  multiplications.  Thus  the  evaluation  of  the  Fourier 

transform  at  all  GH  =  N  points  requires  (G  +  H)N  multiplica- 

2 

tions.  The  conventional  methods  require  N  multiplications. 

If  N  =  G1  G2  G3  G4  ...  G  ,  then  {G±  +  G2  +  Gy  +  G4  +  ...  +  G  )N 

multiples  suffice.  In  particular,  if  N  =  2n,  one  can  com¬ 
pute  the  complete  set  of  Fourier  coefficients  using 
2nN  =  2N  log2  N  multiplications.  For  example,  if  N  =  1024  =  2  ^ , 
then  2N  log2  N  =  20480,  while  =  1048576  which  is  50  times 
20480. 

A2 . 2  Convolution  and  Covariance 

The  convolution  of  two  time  series  in  discrete  form  is 
defined  as  the  inverse  Fourier  transform  of  the  product  of  the 
Fourier  transforms 


N-l 


N-l 


X ( T )  y(t-x)  = 


T  —  0 


N 


X(W)  Y  (W)  e 


2  77  i  t  W 
N 


w=o 


The  convolution  here  must  be  cyclic. 

Sande  (1965)  has  shown  that  the  autocovariance 


R  (x)  = 
xx 


1 

N 


N-t-1 


t=0 


x(t)  x(t+i) 


t  =  0 ,  1 


Lags 


and  the  crosscovariance 


N-t-1 


Rxy  (T)  =  N 


x  ( t )  y(t+x) 


x  =  0 ,  1 


Lags 


t=0 
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can  be  computed  more  economically  by  applying  numerical 
convolution  techniques  rather  than  summing  the  lagged  products. 
But  the  coveriances  are  not  circular  and  the  length  of  the 
series  N  may  not  be  convenient  for  using  Fast  Fourier  transform. 
Both  these  difficulties  can  be  eliminated  by  extending  each  of 
the  series  by  an  appropriate  number  of  zeros.  Thus  extending 
x  series  by  N 1 -N  zeros,  one  gets  a  new  series  such  that 

x'(t)  =  x(t)  t  =  0,  1  ...  N 

x'  (t)  =  0  t  =  N+l  . .  .  N' 


and 

both 


similarly 
series , 

X '  (W) 


for 

using 


y  series.  Then  the 
the  notation  (A2.2) 


N'-l 

SI  x'(t) 
t=0 


,-tw  . 
e  (  jp-  ) 


Fourier  transforms 
can  be  written  as 


of 


Y'  (W) 


N'-l 

y  (s)  e  (  ) 

s=0 


The  inverse  transform  of  the  series  formed  by  the  product  of 
the  Y1  and  the  complex  conjugate  of  X'  is  written  as 


1 

N' 


N  '  1  _^rT 

ET  (X'*  (W)  Y  (W)  )  *  e(  ^-r  )  ■ 

w=o 


1 

N' 


Wx 

N  ' 


C  (t  ) 


N'-l 

El  X'*(W)  Y(W)  e  ( 

w=o 


) 


. 

IX 


i  N-!/'1  n  1  - 1 

N 


=  ff.  21.  ST  H  x-Mt)  y 1  (s )  e(^)  e  <£g)  e(^) 


W=0  t=0  s=0 


1 

N  ' 


N  '  -1  N ' -1  N ' -1 

ZT  JZ  X'*  (t)  y'  (s)  n 
t=0  s=0  W=0 


e  ( 


-W(s-t-x) 

N  ' 


N'-l  N'-l 


-  ^  x'*(t)  y '  (s )  N,  (s-t-x) 


t=0  s=0 


N'-l 


where  N'  6  ,  (s-t-x)  =  2_  e 

t=0 


-2  7t  i  (W-W) 


N' 


where  6N,  (KN)  =  1,  for  integer  K,  otherwise  6  ,  =  0, 


N 

N'-i 

therefore  C(t)  =  XI  x'*(t)  y'(t  +  x). 

t=0 


In  the  time  series  problem  the  series  are  real,  therefore 


N'-l 

C(t)  =  XI  x'(t)  y'(t  +  x). 
t=0 


(A2.5) 


For  x  <N'  -  N,  this  is  just  N  times  the  crosscovariance  R 

xy 

The  same  steps  can  also  be  applied  to  the  autocovariance. 

Thus  the  crosscovariance  (or  autocovariances )  can  be  obtained 
by  taking  the  inverse  Fourier  transform  of  the  products  of 
Fourier  transforms  of  the  series  and  finally  dividing  by  the 
number  of  points  in  the  series. 
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APPENDIX  3 


A3.1.  Apparent  Resistivity  Curves  for  Layered  Anisotropic 

Earth  Models 

In  this  section  the  apparent  resistivity  curves  for 
a  few  3-  and  4-layered  anisotropic  earth  models  are  given 
for  various  values  of  0,  the  angle  from  the  major  axis  of 
anistropy  to  the  axis  of  measurement  x  (Figure  7) . 

The  following  notation  is  used  on  the  plots 
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